/****************************************************************************************************
** Author:  Yasutaka Tominaga 																	   **
** Title:   "Evaluating the Impact of Repeated Leadership Targetings on Militant Group Durability" **
** Journal: International Interactions															   **
** Email:   y-tominaga@yasutakatominaga.com														   **
** Last modified: July 9, 2019																   **
*****************************************************************************************************/

/*For estimating MSM model, I refered to the codes of the presentation slide by Ehsan Karim (2012) (http://fmwww.bc.edu/repec/san2012/karim.san2012.pdf) and Fewell et al. (2004)*/ 
/*Download the replication data, "ReplicationDataSet.dta," to your folder (working directry) */
/*An original Data of leadership targeting comes from Price, B. C. (2012). Targeting top terrorists: How leadership decapitation contributes to counterter- rorism. International Security 36(4), 9–46.*/

/*To open the replication file (.dta), please use the specific version of file matched with your stata version*/
set more off
use "`c(pwd)'/ReplicationDataSet_ver13.dta", clear// for stata version 13

/******************************Variable lists************************************/
/*GroupName: Name of Militant Organizations*/
/*id: group id*/
/*Decap: Leadership Targeting(t)--binary indicator*/
/*Decap_1: previuous value of Decap*/
/*CumDecap: Cumulative Decapitation*/
/*CumDecap_sq: Cumulative Decapitation squared*/
/*count: Numbers of Violent Attacks*/
/*count_0: Baseline value of count*/
/*Intensity: Numbers of Casualties--Killed+Wounded*/
/*Intensity_0: Baseline value of Intensity*/
/*suicide: Numbers of Suicide Attacks*/
/*suicide_0: Baseline value of suicide*/
/*authorities: Numbers of Violent Attacks against authorities*/
/*authorities_0: Baseline value of authorities*/
/*civilians: Numbers of Violent Attacks against targets other than authorities*/
/*civilians_0: Baseline value of civilians*/
/*Ideology_encode: Ideology of each militant groups--category (constant across years)*/
/*Size: Size of each militant group--category (constant across years)*/
/*censor: 1 if observations are censored, 0 otherwise*/
/*EndDummy: 1 if group ends without censoring at specific year, 0 otherwise*/
/*groupend: 1 if an organization ends or not, 0 otherwise*/
/*yrs: Cumulative years*/
/*Start: Start year of militant organizations*/
/*time0: years that observations start*/
/*time1: End year of militant organizations*/
/*f_mil_2: State military capability (Hendrix and Young 2014)*/
/*milcap_0: Baseline value of f_mil_2*/
/*polity2: polity*/
/*polity_0: Baseline value of polity2*/
/*decap_below50th: Cumulative Decapitation before the mid-point of group's lifespan*/
/*decap_over50th: Cumulative Decapitation after the mid-point of group's lifespan*/
/*decap_below50th_sq*/
/*decap_over50th_sq*/


/*****************Calculate Stabilized IPTWs*****************/ 

/*First, calculate the denominator of the equation (1)*/ 
/* Ideology and size are constant across years within groups*/
logistic Decap Decap_1 count Intensity suicide authorities f_mil_2 polity2 i.Ideology_encode i.size
predict pa if e(sample)
replace pa=pa*Decap+(1-pa)*(1-Decap)
sort id time
/*Calculate cumulative probabilities*/
bysort id: replace pa=pa*pa[_n-1] if _n!=1

/*Second, calculate the nominator of the equation (1)*/ 
/*start is the baseline year for each militant group*/
logistic Decap Decap_1 count_0 Intensity_0 suicide_0 authorities_0 milcap_0 polity_0 i.Ideology_encode i.size
predict pa0 if e(sample)
replace pa0=pa0*Decap+(1-pa0)*(1-Decap)
sort id time
bysort id: replace pa0=pa0*pa0[_n-1] if _n!=1

/*Now then we can calcuate SW(t)*/ 
gen w=1/pa
gen sw=pa0/pa


/*****************Calculate Stabilized IPCWs*****************/

/*First, calculate the denominator of the equation (2)*/ 
logistic censor Decap_1 count Intensity suicide authorities i.Ideology_encode i.size polity2 //omit f_mil_2 because f_mil_2 does not have many values at 2008 and so inclusion of it makes it difficult to estimate
predict pcens if e(sample)
replace pcens=pcens*censor+(1-pcens)*(1-censor)
sort id time
bysort id: replace pcens=pcens*pcens[_n-1] if _n!=1 //Calculate cumulative probabilities

/*Second, calculate the nominator of the equation (2)*/ 
logistic censor Decap_1 count_0 Intensity_0 suicide_0 authorities_0 polity_0 i.Ideology_encode i.size
predict pcens0 if e(sample)
replace pcens0=pcens0*censor+(1-pcens0)*(1-censor)
sort id time
bysort id: replace pcens0=pcens0*pcens0[_n-1] if _n!=1 //Calculate cumulative probabilities

/*Now then we can calcuate SW(t) of censoring*/ 
gen csw=pcens0/pcens


/*****************Calculate IPTW*IPCW*****************/ 
gen stabweightcens=sw*csw

//Remove Missing values
foreach v of var * { 
	drop if missing(`v') 
}


/*****************Estimate Causal Impact*****************/ 
/*For each estimate, we get correct standard error through Bootstrap*/

/*Setting Data for Survival Analysis*/
/*allow variable weights within same ID*/
svyset id [pw=stabweightcens]
stset time1 , fail(groupend) origin(Start) enter(time0) exit(time1)

//***************Model 1***************//
/* <The results of Model 1 through Model 3 are shown in Table 2> */
/*Average Treatment Effect conditioning on past history of targeting policies and time-varying confounders*/ 
/*Unweighted*/
streg Decap count_0 Intensity_0 authorities_0 milcap_0 polity_0 i.Ideology_encode i.size, dist(weibull) vce(bootstrap, reps(500) seed(1234))

/*Weighted*/
svy: streg Decap count_0 Intensity_0 authorities_0 milcap_0 polity_0 i.Ideology_encode i.size, dist(weibull)

capture program drop cboot
program define cboot, rclass
svy: streg Decap count_0 Intensity_0 authorities_0 milcap_0 polity_0 i.Ideology_encode i.size, dist(weibull)
return scalar cf = _b[Decap]
end
set seed 123
bootstrap r(cf), reps(500): cboot
estat boot, all


/***************Model 2***************/
/*Cumulative Effect*/
/*Unweighted*/
streg CumDecap count_0 Intensity_0 authorities_0 milcap_0 polity_0 i.Ideology_encode i.size, dist(weibull) vce(bootstrap, reps(500) seed(1234))
/*Weighted*/
svy: streg CumDecap count_0 Intensity_0 authorities_0 milcap_0 polity_0 i.Ideology_encode i.size, dist(weibull)

capture program drop cboot
program define cboot, rclass
svy: streg CumDecap count_0 Intensity_0 authorities_0 milcap_0 polity_0 i.Ideology_encode i.size, dist(weibull)
return scalar cf = _b[CumDecap]
end
set seed 123
bootstrap r(cf), reps(500): cboot
estat boot, all

/***************Model 3***************/
/*Unweighted*/
streg CumDecap CumDecap_sq count_0 Intensity_0 authorities_0 milcap_0 polity_0 i.Ideology_encode i.size, dist(weibull) vce(bootstrap, reps(500) seed(1234))
/*Weighted*/
svy: streg CumDecap CumDecap_sq count_0 Intensity_0 authorities_0 milcap_0 polity_0 i.Ideology_encode i.size, dist(weibull)

capture program drop cboot
program define cboot, rclass
svy: streg CumDecap CumDecap_sq count_0 Intensity_0 authorities_0 milcap_0 polity_0 i.Ideology_encode i.size, dist(weibull)
return scalar cf = _b[CumDecap]
return scalar cf2 = _b[CumDecap_sq]
end
set seed 123
bootstrap r(cf) r(cf2), reps(500): cboot
estat boot, all

/*Marginal Effect*/
/*Left figure in Figure 6*/
svy: streg c.CumDecap##c.CumDecap count_0 Intensity_0 authorities_0 milcap_0 polity_0 i.Ideology_encode i.size, dist(weibull)
predict h0, surv
margins if h0!=., predict(mean) at(CumDecap=(0(0.1)4))
marginsplot, recast(line) recastci(rarea) ciopts(color(gs14)) graphregion(color(white)) plotopts(lcolor(black) lpattern(dash)) title("") ytitle("Predicted Survival Years") xtitle("Number of Leadership Targeting Interventions") saving(f1_a)

/*Show Survival Rate on cumulative effect*/
/*Right figure in Figure 6*/
svy: streg CumDecap CumDecap_sq count_0 Intensity_0 authorities_0 milcap_0 polity_0 i.Ideology_encode i.size, dist(weibull)
stcurve, surv at1(CumDecap=1 CumDecap_sq=1) at2(CumDecap=2 CumDecap_sq=4) at3(CumDecap=3 CumDecap_sq=9) at4(CumDecap=4 CumDecap_sq=16) legend(lab(1 "Targeted=1") lab(2 "Targeted=2") lab(3 "Targeted=3") lab(4 "Targeted=4")) ytitle("Survival Rate") title("") xtitle("Organization Age") graphregion(color(white)) lpattern(solid dot shortdash dash longdash) lcolor(black black black black) saving(f1_b)

/***************Figure 6***************/
gr combine f1_a.gph f1_b.gph, col(2) graphregion(color(white))

/* //For the calculation of minimum/maximum of the quadratic term*/ 
/* svy: streg CumDecap CumDecap_sq count_0 Intensity_0 authorities_0 milcap_0 polity_0 i.Ideology_encode i.size, dist(weibull) nohr*/ 
/* //CumDecap: -.6154466*/ 
/* //CumDecap_sq: .2027314*/ 
/* dis -(-.6154466)/(2*.2027314)*/ 

/***************Model 4***************/
/*  <Table 3> */
/*Cumulative Effec considering when targeting policies had carried out*/ 
/*Cumulative before half years and after half years*/
/*Unweighted*/
streg decap_below50th decap_below50th_sq decap_over50th decap_over50th_sq count_0 Intensity_0 authorities_0 milcap_0 polity_0 i.Ideology_encode i.size, dist(weibull) vce(bootstrap, reps(500) seed(1234))
/*Weighted*/
svy: streg decap_below50th decap_below50th_sq decap_over50th decap_over50th_sq count_0 Intensity_0 authorities_0 milcap_0 polity_0 i.Ideology_encode i.size, dist(weibull)

capture program drop cboot
program define cboot, rclass
svy: streg decap_below50th decap_below50th_sq decap_over50th decap_over50th_sq count_0 Intensity_0 authorities_0 milcap_0 polity_0 i.Ideology_encode i.size, dist(weibull)
return scalar cf = _b[decap_below50th]
return scalar cf2 = _b[decap_below50th_sq]
return scalar cf3 = _b[decap_over50th]
return scalar cf4 = _b[decap_over50th_sq]
end
set seed 123
bootstrap r(cf) r(cf2) r(cf3) r(cf4), reps(500): cboot
estat boot, all

/*For the calculation of minimum/maximum of the quadratic term*/
/* svy: streg decap_below50th decap_below50th_sq decap_over50th decap_over50th_sq count_0 Intensity_0 authorities_0 milcap_0 polity_0 i.Ideology_encode i.size, dist(weibull) nohr*/
/*decap_below50th: -.7053753*/
/*decap_below50th_sq: .3740059*/
/*decap_over50th: .5509746*/
/*decap_over50th_sq: -.6012178*/
/*dis -(-.7053753)/(2*.3740059)*/
/*dis -(.5509746)/(2*-.6012178)*/

/* Figure 7*/
preserve
svy: streg c.decap_below50th##c.decap_below50th c.decap_over50th##c.decap_over50th count_0 Intensity_0 authorities_0 milcap_0 polity_0 i.Ideology_encode i.size, dist(weibull)
predict h2, surv
margins if h2!=., predict(mean) at(decap_below50th=(0(0.1)4))
marginsplot, recast(line) recastci(rarea) ciopts(color(gs14)) graphregion(color(white)) plotopts(lcolor(black) lpattern(dash)) title("") ytitle("Predicted Survival Years") xtitle("Number of Leadership Targeting Interventions" "(Pre-Midpoint of Lifespan)") saving(f2_a)
margins if h2!=., predict(mean) at(decap_over50th=(0(0.1)2))
marginsplot, yscale(range(0 50)) xlabel(0(0.5)2) recast(line) recastci(rarea) ciopts(color(gs14)) graphregion(color(white)) plotopts(lcolor(black) lpattern(dash)) title("") ytitle("Predicted Survival Years") xtitle("Number of Leadership Targeting Interventions" "(Post-Midpoint of Lifespan)") saving(f2_b)
gr combine f2_a.gph f2_b.gph, col(2) graphregion(color(white))
restore

/*****************Testing Covariate Balances between Weighted and Unweighted with the standardized difference*****************/

/***Intensity***/
/*Unweighted*/
reg Intensity Decap Decap_1 count_0 Intensity_0 suicide_0 authorities_0 milcap_0 polity_0 i.Ideology_encode i.size, cluster(id)

/* Calculate CI */
capture program drop cboot
program define cboot, rclass
reg Intensity Decap Decap_1 count_0 Intensity_0 suicide_0 authorities_0 milcap_0 polity_0 i.Ideology_encode i.size, cluster(id)
return scalar cf   = _b[Decap]
return scalar cf2 = _se[Decap]
return scalar cf3 = _b[Decap]/_se[Decap] // the standardized difference
end
set seed 123
bootstrap r(cf) r(cf2) r(cf3), reps(500): cboot
estat boot, all


/*Weighted*/
reg Intensity Decap Decap_1 count_0 Intensity_0 suicide_0 authorities_0 milcap_0 polity_0 i.Ideology_encode i.size [pw=stabweightcens], cluster(id)

/* Calculate CI */
capture program drop cboot
program define cboot, rclass
reg Intensity Decap Decap_1 count_0 Intensity_0 suicide_0 authorities_0 milcap_0 polity_0 i.Ideology_encode i.size [pw=stabweightcens], cluster(id)
return scalar cf   = _b[Decap]
return scalar cf2 = _se[Decap]
return scalar cf3 = _b[Decap]/_se[Decap] // the standardized difference
end
set seed 123
bootstrap r(cf) r(cf2) r(cf3), reps(500): cboot
estat boot, all


/***Authorities***/
/*Unweighted*/
reg authorities Decap Decap_1 count_0 Intensity_0 suicide_0 authorities_0 milcap_0 polity_0 i.Ideology_encode i.size, cluster(id)

/* Calculate CI */
capture program drop cboot
program define cboot, rclass
reg authorities Decap Decap_1 count_0 Intensity_0 suicide_0 authorities_0 milcap_0 polity_0 i.Ideology_encode i.size, cluster(id)
return scalar cf   = _b[Decap]
return scalar cf2 = _se[Decap]
return scalar cf3 = _b[Decap]/_se[Decap] // the standardized difference
end
set seed 123
bootstrap r(cf) r(cf2) r(cf3), reps(500): cboot
estat boot, all

/*Weighted*/
reg authorities Decap Decap_1 count_0 Intensity_0 suicide_0 authorities_0 milcap_0 polity_0 i.Ideology_encode i.size [pw=stabweightcens], cluster(id)

/* Calculate CI */
capture program drop cboot
program define cboot, rclass
reg authorities Decap Decap_1 count_0 Intensity_0 suicide_0 authorities_0 milcap_0 polity_0 i.Ideology_encode i.size [pw=stabweightcens], cluster(id)
return scalar cf   = _b[Decap]
return scalar cf2 = _se[Decap]
return scalar cf3 = _b[Decap]/_se[Decap] // the standardized difference
end
set seed 123
bootstrap r(cf) r(cf2) r(cf3), reps(500): cboot
estat boot, all

/***Suicide***/
/*Unweighted*/
reg suicide Decap Decap_1 count_0 Intensity_0 suicide_0 authorities_0 milcap_0 polity_0 i.Ideology_encode i.size, cluster(id)

/* Calculate CI */
capture program drop cboot
program define cboot, rclass
reg suicide Decap Decap_1 count_0 Intensity_0 suicide_0 authorities_0 milcap_0 polity_0 i.Ideology_encode i.size, cluster(id)
return scalar cf   = _b[Decap]
return scalar cf2 = _se[Decap]
return scalar cf3 = _b[Decap]/_se[Decap] // the standardized difference
end
set seed 123
bootstrap r(cf) r(cf2) r(cf3), reps(500): cboot
estat boot, all

/*Weighted*/
reg suicide Decap Decap_1 count_0 Intensity_0 suicide_0 authorities_0 milcap_0 polity_0 i.Ideology_encode i.size [pw=stabweightcens], cluster(id)

/* Calculate CI */
capture program drop cboot
program define cboot, rclass
reg suicide Decap Decap_1 count_0 Intensity_0 suicide_0 authorities_0 milcap_0 polity_0 i.Ideology_encode i.size [pw=stabweightcens], cluster(id)
return scalar cf   = _b[Decap]
return scalar cf2 = _se[Decap]
return scalar cf3 = _b[Decap]/_se[Decap] // the standardized difference
end
set seed 123
bootstrap r(cf) r(cf2) r(cf3), reps(500): cboot
estat boot, all


/***Count***/
/*Unweighted*/
reg count Decap Decap_1 count_0 Intensity_0 suicide_0 authorities_0 milcap_0 polity_0 i.Ideology_encode i.size, cluster(id)

capture program drop cboot
program define cboot, rclass
reg count Decap Decap_1 count_0 Intensity_0 suicide_0 authorities_0 milcap_0 polity_0 i.Ideology_encode i.size, cluster(id)
return scalar cf   = _b[Decap]
return scalar cf2 = _se[Decap]
return scalar cf3 = _b[Decap]/_se[Decap] // the standardized difference
end
set seed 123
bootstrap r(cf) r(cf2) r(cf3), reps(500): cboot
estat boot, all

/*Weighted*/
reg count Decap Decap_1 count_0 Intensity_0 suicide_0 authorities_0 milcap_0 polity_0 i.Ideology_encode i.size [pw=stabweightcens], cluster(id)

capture program drop cboot
program define cboot, rclass
reg count Decap Decap_1 count_0 Intensity_0 suicide_0 authorities_0 milcap_0 polity_0 i.Ideology_encode i.size [pw=stabweightcens], cluster(id)
return scalar cf   = _b[Decap]
return scalar cf2 = _se[Decap]
return scalar cf3 = _b[Decap]/_se[Decap] // the standardized difference
end
set seed 123
bootstrap r(cf) r(cf2) r(cf3), reps(500): cboot
estat boot, all

/*Figure 5*/
/*Values above stored in Figure3.dta file separately and re-imported into Stata for Plot*/ 
use "`c(pwd)'/Figure3Data_ver13.dta", replace
graph box value, over(var) over(weight) asyvars intensity(80) ytitle("Standardized Coefficients") scheme(s2mono) graphregion(color(white))


/*****************Other Figures*****************/ 
/*Figure 1*/
use "`c(pwd)'/Figure1Data_ver13.dta", replace
graph hbar Strikes, bargap(0) over(Group_01, label(labsize(vsmall))) over(Names, sort(Strikes) descending label(labsize(vsmall))) asyvars bar(1, color(black*1)) bar(2, color(grey*0.65)) bar(3, color(grey*0.25)) bar(4, color(grey*0.1)) graphregion(color(white) margin(small)) intensity(70)  legend(position(5) ring(0) size(1)) nofill ytitle("Numbers of Drone Strikes")

/*Figure 2*/
use "`c(pwd)'/Figure2Data_ver13.dta", replace
graph bar SeniorExpertranking others, over(Targetsorganization, sort(1) label(labsize(vsmall)))  graphregion(color(white) margin(small)) intensity(70) stack bar(1, color(black*0.8)) bar(2, color(black*0.4)) ytitle("Number of Targetings") legend(label(1 "Senior level militants") label(2 "Others"))

/*Figure 3 and Figure 4 are created using Adobe Illusrator*/
/*Figure 8 is created using ArcGIS 10.4.1*/


//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////Online Appendix///////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////


////////////////////////////Include civilian target as a time-varying confounders rather than violent attacks against authorities////////////////////////////
//**************************Repeat the same process using civilians_0 and civilians***********************************************//
use "`c(pwd)'/ReplicationDataSet_ver13.dta", replace

/*First, calculate the denominator of the equation (1)*/ 
logistic Decap Decap_1 count Intensity suicide civilians f_mil_2 polity2 i.Ideology_encode i.size // Ideology and size are constant across years within groups
predict pa if e(sample)
replace pa=pa*Decap+(1-pa)*(1-Decap)
sort id time
bysort id: replace pa=pa*pa[_n-1] if _n!=1 //Calculate cumulative probabilities

/*Second, calculate the nominator of the equation (1)*/ 
logistic Decap Decap_1 count_0 Intensity_0 suicide_0 civilians_0 milcap_0 polity_0 i.Ideology_encode i.size //start is the baseline year for each militant group
predict pa0 if e(sample)
replace pa0=pa0*Decap+(1-pa0)*(1-Decap)
sort id time
bysort id: replace pa0=pa0*pa0[_n-1] if _n!=1

/*Now then we can calcuate SW(t)*/ 
gen w=1/pa
gen sw=pa0/pa


/*****************Calculate Stabilized IPCWs*****************/

/*First, calculate the denominator of the equation (2)*/ 
logistic censor Decap_1 count Intensity suicide civilians i.Ideology_encode i.size polity2 //omit f_mil_2 because f_mil_2 does not have many values at 2008 and so inclusion of it makes it difficult to estimate
predict pcens if e(sample)
replace pcens=pcens*censor+(1-pcens)*(1-censor)
sort id time
bysort id: replace pcens=pcens*pcens[_n-1] if _n!=1 //Calculate cumulative probabilities

/*Second, calculate the nominator of the equation (2)*/ 
logistic censor Decap_1 count_0 Intensity_0 suicide_0 civilians_0 polity_0 i.Ideology_encode i.size
predict pcens0 if e(sample)
replace pcens0=pcens0*censor+(1-pcens0)*(1-censor)
sort id time
bysort id: replace pcens0=pcens0*pcens0[_n-1] if _n!=1 //Calculate cumulative probabilities

/*Now then we can calcuate SW(t) of censoring*/ 
gen csw=pcens0/pcens


/*****************Calculate IPTW*IPCW*****************/ 
gen stabweightcens=sw*csw

//Remove Missing values
foreach v of var * { 
	drop if missing(`v') 
}

/*****************Estimate Causal Impact*****************/ 
/*For each estimate, we get correct standard errors through Bootstrap*/

/*Setting Data for Survival Analysis*/
svyset id [pw = stabweightcens] //allow variable weights within same ID
stset time1 , fail(groupend) origin(Start) enter(time0) exit(time1)

/*Model 1*/
/*<<Table 3 in Appendix from Model 1 through Model 3>>*/
/****Average Treatment Effect conditioning on past history of targeting policies and time-varying confounders*** */ 
/*Unweighted*/
streg Decap count_0 Intensity_0 civilians_0 milcap_0 polity_0 i.Ideology_encode i.size, dist(weibull) vce(bootstrap, reps(500) seed(1234))

/*Weighted*/
svy: streg Decap count_0 Intensity_0 civilians_0 milcap_0 polity_0 i.Ideology_encode i.size, dist(weibull)

capture program drop cboot
program define cboot, rclass
svy: streg Decap count_0 Intensity_0 civilians_0 milcap_0 polity_0 i.Ideology_encode i.size, dist(weibull)
return scalar cf = _b[Decap]
end
set seed 123
bootstrap r(cf), reps(500): cboot
estat boot, all

/*Cumulative Effect*/
/*Model 2*/
/*Unweighted*/
streg CumDecap count_0 Intensity_0 civilians_0 milcap_0 polity_0 i.Ideology_encode i.size, dist(weibull) vce(bootstrap, reps(500) seed(1234))
/*Weighted*/
svy: streg CumDecap count_0 Intensity_0 civilians_0 milcap_0 polity_0 i.Ideology_encode i.size, dist(weibull)

capture program drop cboot
program define cboot, rclass
svy: streg CumDecap count_0 Intensity_0 civilians_0 milcap_0 polity_0 i.Ideology_encode i.size, dist(weibull)
return scalar cf = _b[CumDecap]
end
set seed 123
bootstrap r(cf), reps(500): cboot
estat boot, all

/*Model 3*/
/*Unweighted*/
streg CumDecap CumDecap_sq count_0 Intensity_0 civilians_0 milcap_0 polity_0 i.Ideology_encode i.size, dist(weibull) vce(bootstrap, reps(500) seed(1234))
/*Weighted*/
svy: streg CumDecap CumDecap_sq count_0 Intensity_0 civilians_0 milcap_0 polity_0 i.Ideology_encode i.size, dist(weibull)

capture program drop cboot
program define cboot, rclass
svy: streg CumDecap CumDecap_sq count_0 Intensity_0 civilians_0 milcap_0 polity_0 i.Ideology_encode i.size, dist(weibull)
return scalar cf = _b[CumDecap]
return scalar cf2 = _b[CumDecap_sq]
end
set seed 123
bootstrap r(cf) r(cf2), reps(500): cboot
estat boot, all

/*Model 4*/
/*Table 4 in Online Appendix*/
/*Cumulative Effec considering when targeting policies had carried out*/ 
/*Cumulative before half years and after half years*/
/*Unweighted*/
streg decap_below50th decap_below50th_sq decap_over50th decap_over50th_sq count_0 Intensity_0 civilians_0 milcap_0 polity_0 i.Ideology_encode i.size, dist(weibull) vce(bootstrap, reps(500) seed(1234))
/*Weighted*/
svy: streg decap_below50th decap_below50th_sq decap_over50th decap_over50th_sq count_0 Intensity_0 civilians_0 milcap_0 polity_0 i.Ideology_encode i.size, dist(weibull)

capture program drop cboot
program define cboot, rclass
svy: streg decap_below50th decap_below50th_sq decap_over50th decap_over50th_sq count_0 Intensity_0 civilians_0 milcap_0 polity_0 i.Ideology_encode i.size, dist(weibull)
return scalar cf = _b[decap_below50th]
return scalar cf2 = _b[decap_below50th_sq]
return scalar cf3 = _b[decap_over50th]
return scalar cf4 = _b[decap_over50th_sq]
end
set seed 123
bootstrap r(cf) r(cf2) r(cf3) r(cf4), reps(500): cboot
estat boot, all


/**************************************TABLE 4 in the Main Article **************************************/

/*** Weighted ***/
svy: streg decap_below50th c.decap_below50th##c.decap_below50th decap_over50th c.decap_over50th##c.decap_over50th count_0 Intensity_0 authorities_0 milcap_0 polity_0 i.Ideology_encode i.size, dist(weibull) nohr
predict h2, surv

/*** Al-Qaeda ***/
margins if h2!=., predict(mean) at(decap_below50th=0 decap_over50th=3 count_0=0 Intensity_0=0 authorities_0=0 milcap_0=0.777 polity_0=8 Ideology_encode=4 size=1000)

/*** Haqqani Network ***/
margins if h2!=., predict(mean) at(decap_below50th=0 decap_over50th=1 count_0=0 Intensity_0=0 authorities_0=0 milcap_0=0.602 polity_0=-7 Ideology_encode=4 size=1000)

/*** IMU ***/
margins if h2!=., predict(mean) at(decap_below50th=0 decap_over50th=1 count_0=0 Intensity_0=0 authorities_0=0 milcap_0=0.993 polity_0=7 Ideology_encode=4 size=1000)

/*** Taliban ***/ 
margins if h2!=., predict(mean) at(decap_below50th=0 decap_over50th=1 count_0=0 Intensity_0=0 authorities_0=0 milcap_0=0.941 polity_0=8 Ideology_encode=4 size=1000)

/*** UnWeighted ***/
streg decap_below50th c.decap_below50th##c.decap_below50th decap_over50th c.decap_over50th##c.decap_over50th count_0 Intensity_0 authorities_0 milcap_0 polity_0 i.Ideology_encode i.size, dist(weibull) nohr
predict h3, surv

/*** Al-Qaeda ***/
margins if h3!=., predict(mean) at(decap_below50th=0 decap_over50th=3 count_0=0 Intensity_0=0 authorities_0=0 milcap_0=0.777 polity_0=8 Ideology_encode=4 size=1000)

/*** Haqqani Network ***/
margins if h3!=., predict(mean) at(decap_below50th=0 decap_over50th=1 count_0=0 Intensity_0=0 authorities_0=0 milcap_0=0.602 polity_0=-7 Ideology_encode=4 size=1000)

//*** IMU ***/
margins if h3!=., predict(mean) at(decap_below50th=0 decap_over50th=1 count_0=0 Intensity_0=0 authorities_0=0 milcap_0=0.993 polity_0=7 Ideology_encode=4 size=1000)

/*** Taliban ***/ 
margins if h3!=., predict(mean) at(decap_below50th=0 decap_over50th=1 count_0=0 Intensity_0=0 authorities_0=0 milcap_0=0.941 polity_0=8 Ideology_encode=4 size=1000)


/*****************Robustness Check*****************/ 
/* Use Marginal Structural Cox Hazard Model */

/*** Table 1 in Online Appendix (Model 1 through Model 3) ***/
/*** Re-estimate Model 1 ***/
/*  Unweighted */ 
stcox Decap count_0 Intensity_0 authorities_0 milcap_0 polity_0 i.Ideology_encode i.size, vce(bootstrap, reps(500) seed(1234))
/*  Weighted */
svy: stcox Decap count_0 Intensity_0 authorities_0 milcap_0 polity_0 i.Ideology_encode i.size

capture program drop cboot
program define cboot, rclass
svy: stcox Decap count_0 Intensity_0 authorities_0 milcap_0 polity_0 i.Ideology_encode i.size
return scalar cf = _b[Decap]
end
set seed 123
bootstrap r(cf), reps(500): cboot
estat boot, all

/** Re-estimate Model 2 **/
/*  Cumulative Effect  */
/*  Unweighted  */
stcox CumDecap count_0 Intensity_0 authorities_0 milcap_0 polity_0 i.Ideology_encode i.size, vce(bootstrap, reps(500) seed(1234))
/*  Weighted  */
svy: stcox CumDecap count_0 Intensity_0 authorities_0 milcap_0 polity_0 i.Ideology_encode i.size

capture program drop cboot
program define cboot, rclass
svy: stcox CumDecap count_0 Intensity_0 authorities_0 milcap_0 polity_0 i.Ideology_encode i.size
return scalar cf = _b[CumDecap]
end
set seed 123
bootstrap r(cf), reps(500): cboot
estat boot, all

/*** Re-estimate Model 3 ***/
/*  Including Squared Term */
/*  Unweighted  */
stcox CumDecap CumDecap_sq count_0 Intensity_0 authorities_0 milcap_0 polity_0 i.Ideology_encode i.size, vce(bootstrap, reps(500) seed(1234))
/*  Weighted  */
svy: stcox CumDecap CumDecap_sq count_0 Intensity_0 authorities_0 milcap_0 polity_0 i.Ideology_encode i.size

capture program drop cboot
program define cboot, rclass
svy: stcox CumDecap CumDecap_sq count_0 Intensity_0 authorities_0 milcap_0 polity_0 i.Ideology_encode i.size
return scalar cf = _b[CumDecap]
return scalar cf2 = _b[CumDecap_sq]
end
set seed 123
bootstrap r(cf) r(cf2), reps(500): cboot
estat boot, all

/*** Re-estimate Model 4 ***/
/*** Table 2 in Online Appendix ***/
/*  Cumulative before half years and after half years  */
/*  Unweighted  */
stcox decap_below50th decap_below50th_sq decap_over50th decap_over50th_sq count_0 Intensity_0 authorities_0 milcap_0 polity_0 i.Ideology_encode i.size, vce(bootstrap, reps(500) seed(1234))
/*  Weighted  */
svy: stcox  decap_below50th decap_below50th_sq decap_over50th decap_over50th_sq count_0 Intensity_0 authorities_0 milcap_0 polity_0 i.Ideology_encode i.size

capture program drop cboot
program define cboot, rclass
svy: stcox  decap_below50th decap_below50th_sq decap_over50th decap_over50th_sq count_0 Intensity_0 authorities_0 milcap_0 polity_0 i.Ideology_encode i.size
return scalar cf = _b[decap_below50th]
return scalar cf2 = _b[decap_below50th_sq]
return scalar cf3 = _b[decap_over50th]
return scalar cf4 = _b[decap_over50th_sq]
end
set seed 123
bootstrap r(cf) r(cf2) r(cf3) r(cf4), reps(500): cboot
estat boot, all


///////////////////////////////////////////////////////////////////////////////END//////////////////////////////////////////////////////////////////////////////////////
